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1.  INTRODUCTION 

This  report  describes  the  implementation  of  a  Third  order  Butterworth  filter  in  the 
Advanced  Continuous  Simulation  Language  (ACSL)  [1].  Transient  responses  of  the  filter  in 
the  time  domain,  where  the  time  domain  expressions  were  derived  from  both  s  and  i  domain 
representations,  are  observed  and  compared  with  the  responses  from  the  analytical  time 
domain  expression  for  the  filter.  Inputs  are  step  and  sine  wave  functions.  There  are  two  z 
domain  implementations  used:  the  bilinear  transformation  method  and  the  impulse  invariant 
method  [2].  These  are  both  based  on  transforming  the  continuous  representation  of  the  filter 
in  the  s  domain  to  the  z  domain. 

An  ACSL  macro  for  implementing  Z*transfer  functions  was  developed  and  used  in  the 
transient  response  testing.  The  macro  algorithm  is  based  on  the  “M-Method”  canonical  form 
of  the  transfer  function.  Only  zero  initial  conditions  are  allowed  for  in  the  algorithm. 

Bode  plots  of  the  s  and  Z  transfer  functions  were  made  with  an  ACSL  program.  The 
program  and  Bode  plots  are  both  described  in  this  report. 

2.  MATH  MODELS  FOR  A  THIRD  ORDER  BUTTERWORTH  FILTER 
IN  THE  Z  AND  s  DOMAINS 

A  third  order  Butterworth  filter  in  its  normalized  s  domain  form  looks  as  follows 


s3  +  2s2  +  2s  +  1 


or, 


_ 1 _ 

(s  +  1)  (s2  +  S  +  1) 

If  the  cutoff  frequency  a*  is  shown,  the  filter  takes  the  form 


(-  +  1) 
wo 


(2) 


The  s  domain  representation  of  the  output  of  the  filter  for  a  step  input  of  unity  gain  is 


a) 


2  2 
S  (S  +  WQ)  (S  +  WqS  +  WQ) 


This  can  also  be  written  as 
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which  is  easily  transferred  to  its  equivalent  in  the  time  domain, 

u(t)  -  e"ta)ot  -  ^  V^sin  (  uQt) 


(3) 


where  u(t)  is  a  step  function  with  unity  gain.  Similarly  for  a  unity  gain  sine  wave  of  frequency 
yr  the  time  domain  output  of  the  filter  is 
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These  analytical  time  domain  expressions  will  be  used  to  compare  with  the  s  and  Z  domain 
results  in  the  transient  response  simulation. 

The  s  domain  filter  expression  used  in  the  simulation  is,  from  (2), 


3 - 5 - 5  3 

8  +  2(1)  8  +  2d)  S  +  (•) 

o  o  o 


In  the  impulse  invariant  method  of  transforming  a  transfer  function  from  the  s  domain  to  the 
Z  domain,  the  impulse  response  of  the  transfer  function  is  mapped  directly  into  the  Z  domain. 
Equation  (5)  can  be  rewritten  as 
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<D  8  +  0. 5  (D 
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0.5  u) 


w  (R0  +  R^z”1  +  R2z-2) 


-1  -2  -3 

So  +  +  S2z  +  S3z 


where 


Z  =  e 


T  is  the  sampling  interval 


or 

"  r  T  /3 

u  =  e  cos  (  a)  T) 

2  o 


3T~  1  Si 

v  =  e  sin  (  <jj_T) 


2.2  -oT 
u  +  v  =  e 
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R  =  -U  +  —  +  U2  +  V2 

1  /3 


R-  =  (u2  +  V2)  (1  -  U  -  ) 

^  /3 


So  =  1 


S.^  =  -  (2u  +  u2  +  v2) 


S2  =  (u2+v2) (l+2u) 


Upon  implementing  (6)  one  immediately  notices  the  gain  is  <u0/T  rather  than  1.  This  is  due 
in  part  to  the  fact  that 

00 

X(z)  =  -TJ-  I  X  (s  +  j  k) 

k=-°° 

2 

where  X  is  an  arbitrary  transfer  function  [3]  and  SK.  x  (s+j  K)  is  the  repeated  spectrum 
(due  to  the  sampling  of  X(t))of  the  continuous  s-transfer  function.  Since  the  continuous 
representation  of  X  in  the  s  domain  repeats,  one  must  use  care  in  selecting  T.  Hopefully  the 
frequency  content  of  the  input  signal  is  known  to  aid  in  this  determination. 

In  the  bilinear  transformation  method  s  is  directly  replaced  by  z-l/z+1  in  Equation  (5). 
This  change  of  variable  maps  the  imaginary  axis  of  the  s-plane  to  the  unit  circle  of  the  7-plane 
[21  eq,is  replaced  by  tan  <uo  T/  2  in  (5).  The  analogy  in  the  frequency  domain  being  s— ja>,  then 
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Thus  Equation  (S)  becom. 


aQ  =  a3/B 
al  =  a2  =  3ao 


bQ  *3 

bx  -  -  (3+2a  -  2a2  -  3a3)/0 

b2  -  +  (3-2a  -  2a2  +  3a3)/8 

b3  -  -  (l-2a  +  2a2  -  a3)/g 

CD T 

a  *  ban  — j — 

B  =  1  +  2a  +  2a2  +  a3 

2ir 

Although  this  transfer  function  representation  also  repeats  itself  every  “j~K  rad  /  sec,  there  is 
no  degradation  in  the  filter  roll-off  characteristics  as  will  be  shown  in  Section  3. 

As  in  the  s  domain  there  are  canonical  forms  for  Z-transfer  functions.  The  “M-Method” 
canonical  implementation  of  the  third  order  Butterworth  filter  is  shown  in  Figure  1.  The  figure 
illustrates  the  representation  for  Equation  (7).  Equation  (7)  becomes 

M=  ( INPUT  -  b1Z_1M  -  b2Z"2M  -  b3Z*3M)  (8) 

OUTPUT  -  (aQ  +  a^Z-1  +  a2z”2  +  a3Z_3)  M  (9) 

ACSL  macro  ZXSFRM  was  designed  to  implement  a  generalized  form  of  the  M-method  that 
can  be  used  for  any  order  Z-transform  whose  highest  power  in  the  numerator  does  not  exceed 
that  of  the  denominator.  The  macro  algorithm  is  described  in  Appendix  A 


12 


The  continuous  •  domain  filter  representation  was  also  implemented  in  the  time  domain 
through  the  RK2  numerical  integration  scheme  with  1/s  interpreted  as  an  integration 
operation.  The  ACSL  macro  TRANS  [1]  was  used  to  implement  the  numerical  integration  of 
the  filter  states.  TRANS  contains  an  algorithm  which  generates  the  state  equations  based  on 
the  M -method  canonical  form  shown  in  Figure  2  for  Equation  (S).  It  should  be  noted  that  non¬ 
zero  intitial  conditions  on  the  states  of  the  filter  are  not  available  in  TRANS  or  ZXSFRM. 

The  direct  numerical  integration  of  the  s  domain  filter  u  11  be  compared  with  the  analytical 
response  and  the  two  Z  domain  implementations. 


UTPUT 
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3.  BODE  PT  "T  PROGRAM  AND  RESULTS  FOR  A  THIRD  ORDER 
BUTTF  -  <  ”*TH  FILTER  IN  THE  Z  AND  i  DOMAINS 

In  order  tc  .ec  numerical  oosfflrients  of  the  impulse  invariant  ana  bilinear 
tranefonnation  tw*r.  of  the  digital  Putter* orth  filter,  an  ACSL  Bode  plot  program  was 
«o  produce  z«traniform  eloti.  The  orldiil  alnarithm  utttlaeB  an  wwbiHtd 
representation  of  the  polynomials  which  is  weB  soiled  for  DO  loop  Implementation.  Figure  3 
ie  a  listing  of  the  Bode  plot  program  whl  throe  macros  added  to  it  for  generation  of  the 
ooeMelcn  of  the  three  filters;  the  i  do— in  vomUn.CCN  IB,  and  the  twoZdomatn  versions, 
CWVEMT  and  CTILNR.  The  macros  are  weed  to  eelenlme  the  cotffiri— s  for  different  break 
frequencies  and  sampling  Intervals  In  a  farm  that  would  not  stutter  ap  what  attempts  to  be  a 
g—ral  purpose  Bode  plot  program.  The  aquarians  in  CNVRNT,  CB1LNR  and  CCNTS  ate 
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PROGRAM  BODE  PLOTS  FOR  TPAnSFER  FNS  IN  THE  S  AnO  Z  DOMAINS 
••  DAVID  B.  MERRIMAN 

••  12  OCT  78 

"  REDSTONE  ARSENAL*  AL 

"  BASED  ON  A  PROGRAM  WRITTEN  BY 

••  E.E.L.  MITCHELL  AND  AN  ALGORITHM  GIVEN  BY 

"  MERV  BUDGE. 

MACRO  CNVRNT<U*V.WB»OT*KG) 

MACRO  REDEFINE  K.E.CtS 
PROCEDURAL  <U*V»WB«OT.KG> 

E*EXP(-WB«OT> 

C=SQRT  <E)*COS(SQRT (3. ) •O.S*W0*OT) 

S=SORT  <E) •SINISORT (3. ) »0.5*WB*OT) /SORT « 3. > 

U(3)*0.0 

UC2»*WB«(-C*S*E> 

un>*WB*E«<I.-C-S) 

V<4)*1.0 

V<3)*-2.«C-E 

V<2)®*E*<1.*2<*C» 

V(l)a-E**2 

K»(U(1)*U(2))/(V(1)«V(2>«V(3>*V<4>> 

U ( I ) «U I 1 » /K 
U  <  2 ) *U ( 2 ) /K 
U(3)*U(3)/K 
FNO 

MACRO  END 

MACRO  CRILNR<R.S«WB.OT*MG) 

MACRO  REDEFINE  AL.OEN 
PROCEDURAL (R . S*WB  «OT  «  KG ) 

AL=TAN(O.S*WB«OT) 

DEN=1.*2.*AL*2.*AL»*24AL»«3 

R<4)*AL**3/DEN 

P(3>»3.«AL**3/DEN 

R<2)*P<3) 

P ( 1 ) *R (4) 

S (4) »1 .0 

S<3)»<-3.-2.*AL*2.<*AL**2*3.«AL»*3>/DEN 
S12)«(3.-2.*AL-2.*AL»»20.*AL»*3>/DEN 
S ( 1 ) « < -I . ♦2.*AL-2.*AL*»2*AL«t3) /DEN 
END 

MACRO  ENO 

MACRO  CCNTS(A.B.WB) 

A(I)*WB*»3 
Btl)*l.O 
B(2)*2.0«WB 
B<3)*»2.0*WB**2 
B  <4 ) *A ( I ) 

MACRO  ENO 


••LIST  OF  SYM80LS" 


"A-NUMERATOR  COEFFICIENT  ARRAY  IN  DECREASING  POWERS  OF  S  OR  Z**(-l>« 
"B-OENOMINATOR  COEFFICIENT  ARRAY  IN  DECREASING  POWERS  OF  S  OR  Z** (-!>•• 
••CNTNUS-IF  .T,  XSFER  FN  IN  S»  IF  ,F.  XSFER  FN  IN  Z“ 

••DEN-DENONINATOR  OF  G« 

••G-TRANSFER  FN  W/A  CONST  GAIN  (SEE  KG>" 

••OT-SAMPLING  PERIOD  FOR  Z-TRANSFORM** 

••GAIN-GAIN  OF  THE  TRANSFER  FUNCTION  (DB>»* 

••KW-GEOMETRIC  PROGRESSION  MULTIPLYING  FACTOR,  W(I,l)sW(I)*KW** 
••KG-OPTIONAL  INPUT  FOR  THE  CONST  GAIN  OF  THE  XSFER  FN" 

••M-ORDER  OF  THE  NUMERATOR  POLYNOMIAL  ♦  !•• 

••N-ORDER  OF  THE  DENOMINATOR  POLYNOMIAL  ♦  1'* 

••NO-TOTAL  NUMBER  OF  FREQUENCY  POINTS  calculated  BETWEEN  WMN  AND  WMX" 
'•NUN-NUMERATOR  of  g»* 

••PHASE-PHASE  OF  THE  TRANSFER  FUNCTION  (DEG)** 

••LOGW-ALOGIO(W)** 

'•MAG-CALCULATED  GAIN  of  THE  XSFER  FN** 

••W-PRESENT  VALUE  OF  FREQUENCY  (RAO/SEC)** 

"WCPS-PRESENT  VALUE  OF  FREQUENCY  (CPS)1* 

••WMN-LOWEST  FREQUENCY  (RAO/SEC)*' 

••WMX-HIGHEST  FREQUENCY  (RAD/SEC) •< 

"X-REAL  PART  OF  G" 

"Y-IMAGINARY  PART  OF  G“ 

••INPUTS  M,  N,  WMN,  WMX,  MPTS,  A,  B,  DT »  KG** 

••OUTPUTS  LOGW,  W,  WCPS*  GAIN,  PHASE,  MAG,  X,  y" 


•• - — — - - - ———GLOBAL  CONSTANTS 

LOGICAL  OUMP 

CONSTANT  RMN  «  1.0E-30,  RMX  *  I.OE30  ,  RADDEG.57, *957795 
,  PI*3.1<*159?65*  DUMP  *  .FALSE. 

TWOPI  «  Z.O*P I 
INITIAL 

REAL  A (SO)  ,  B (50) 

INTEGER  M  ,  N 

LOGICAL  CNTNUS 

CONSTANT  M  ■  1  ,  N  >  1  ,  CNTNUS  *  .TRUE. 

,  A  >  50*0.0  ,  B  ■  50*0.0 

CONSTANT  OT  »  0*010 

CONSTANT  KG  >  1.0 

•• - - - - - — - - - FREQUENCY  IS  SWEPT  IN  A  GEOMETRIC 

PROGRESSION  OVER  A  RANGE  FROM  *w  MN*  TO  *W  MX*,  *MPTS*  POINTS 

PLOT  ANO  KW  IS  THE  MULTIPLIER  BETWEEN  FREQUENCIES - 

CONSTANT  WMN  >1.0  ,  WMX  >  350,  ,  MPTS  =  50.0 

KW  *  10.0** (ALOG10 (WMX/WMN) /MPTS) 

W  ■  WMN 


FlguroS.  (ConMnutd) 
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" - ..... - - .^-SET  BREAK  FREQUENCY  (RAO)  AND  CALCULATE..* 

XSFER  fn  COEFFICIENTS  ■ 

CONSTANT  WBCpS  ■  7*0 
MB  «  WBCPS* TWORI 

LOGICAL  BILNG  ,  NyRNT 

CONSTANT  BILNR**FAL$E*t  NyRNT-.FALSE. 

IF(.N0Ti|BILNR)G0  TO  SN099 
CBILNRIA*  B"WB*  OT*  KB) 

SN999.*IF(.N0T.NVRNT>G0  TO  SN998 
CNVRNT ( A  *  B«NB*  OT*  KG) 

SN998* .CONTINUE 

IF(.NOT%CNTNUS)GO  TO  SN997 
CCNTSM*  B«WB) 

SN997*. CONTINUE 
END  S  "INITIAL" 

DYNAMIC 

- - - ......... - CALCULATE  REAL  A  NO  IMAGINARY  PARTS  OF  G  " 

IF (CNTNUS)CALL  BODES 
IFf.NOT\)CNTNUS)CALL  BOBEZ 

- - - - CALCULATE  GAIN  AND  PRASE  OF  G  " 

MAG  ■  SORT (X**2  ♦  y«*2> 

GAIN  •  20**ALOG10 (NAG) 

PHASE  >  RAOOEG*ATRN2(Y*X) 

.............. ——USE  LOGARITHM  OF  FREQUENCY  FOR  PLOTTING  " 

LOOM 


ALOGli(W) 


W/TWOPI 


V  «  KW*W 

TERMT ( Wf 6T *MMX I 
ENO  S  "OYNAMIC" 

TERMINAL 

IF (DUMP) CALL  DEBUG 
END  S  "TERMINAL" 

END  $  "PROGRAM* 


-MAKE  FREQUENCY  AVAILtBLE  IN  HERTZ 
-ADVANCE  FREQUENCY  GEOMETRICALLY 


direct  implementations  of  equations  (6),  (7)  and  (5)  of  Section  2,  respectively.  hitheCNVRNT 
macro,  instead  of  modifying  the  gain  of  the  impulse  invariant  Z-transform  by  multiplying  by 
T,  the  sampling  interval  (see  Section  2),  the  gain  of  the  transform  at  <u=0 
was  used: 


For  s=j<w 


Z 


e'jwT 


and  at  <o=0 

Z  =*  1 


Hence  for  a  generalized  Z  transform 

Output  =  ao  *  alZ  1  *  a2Z  2 

Input  ”  +  ^  jj-1  +  ^  z“2  + 

r>  1  u2 

the  gain  at  tm= 0  is 

Output  =  ao  *  al  +  a2  + 

Input  bQ  +  +  b2  +  ... 


From  Equation  (6) 


Output 

Input 


3.161xl03  »  gain 


(10) 


when 


u  *  0 

<a  *  5HZ  =  lOtr  rad/sec 
o 

T  *  0.01  sec 

\  (Dg 

Note  that  the  gain  is  not  “p  instead  the  gain  is  equivalent  to  the  break  frequency  (rad) 
divided  by  the  sampling  interval.  Perhaps  some  insight  may  be  gained  as  to  the  origin  of  the  cue 


18 


term  by  observing  that  the  impulse  response  of  the  continuous  filter  in  the  time  domain 
also  contains  a  gain  factor  of  coo.  From  Equation  (5),  the  impulse  response  of  the  continuous 
Butterworth  filter  in  the  time  domain  is, 


fvl 


-  e 


rjf 

cos 


<l) 


/r 


U>r>t  + 


/r 


sin 


(11) 


Next  in  Figure  3  is  a  definition  of  the  standard  program  variables  along  with  a  listing  of  the 
input  and  output  variables.  In  the  prelNITIAL  section,  global  constants  are  defined.  These 
are  parameters  such  as  computer  arithmetic  minimum  and  maximum,  RMN  and  RMX;  rr;  the 
conversion  factor  for  radians  to  degrees,  RADDEG;  and  the  logical  variable  DUMP  which  is 
used  to  control  the  debug  printing  of  the  program  variables  for  the  final  communication 
interval. 


In  the  INITIAL  section  the  parameters  for  this  particular  program  are  defined.  For  the 
Bode  plots  the  frequency  w  is  swept  in  a  geometric  progression.  The  plotted  points  will  then  be 
equally  spaced  when  the  abscissa  is  logiow.  The  multiplier  between  frequencies,  kW,  is 
calculated  based  on  the  minimum  and  maximum  values  for  <u,  WMN  and  WMX,  and  the 
number  of  plotted  points  desired,  MPTS.  w  is  initialized  to  WMN. 

The  code  following  the  determination  of  KW  and  ending  with  the  INITIAL  section  is 
specially  added  for  the  Butterworth  filter  plots.  The  break  frequency  in  Hz,  WBCPS,  is  set  by  a 
CONSTANT  statement  and  then  converted  to  radians  as  represented  by  the  variable  WB. 
BILNR  and  NVRNT  are  defined  to  be  logical  variables  which,  when  TRUE,  cause  the 
calculation  of  the  bilinear  Z-transform  and  the  impulse  invariant  transform  coefficients, 
respectively.  If  CNTNUS  is  TRUE  the  coefficients  for  the  s  domain  Butterworth  filter  are 
calculated.  6bviously  only  one  of  three  logical  variables  CNTNUS,  BILNR  and  NVRNT 
should  be  TRUE  for  a  particular  run. 

The  DYNAMIC  section  contains  the  code  for  calculating  the  magnitude  and  phase  of  the 
transfer  function.  The  real  and  imaginary  components  of  the  transfer  function,  x  and  y,  are 
calculated  in  FORTRAN  subroutine  BODES  which  has  an  entry  BODEZ  for  Z-transfcr 
functions.  For  the  abscissa  the  log  of  w  is  used, 


LOGW  -  ALOGIO(W) 


(12) 


a)  is  made  available  in  Hz  for  printer  output, 

WCPS=W/TWOPI. 

to  is  incremented  and  the  DYNAMIC  section  is  re-executed  until  <o>WMX, 

W=KW*W 

TERMT  (W.GT.WMX) 

There  are  no  integrations  made  so  a  DERIVATIVE  section  is  not  needed. 

Figure  4  is  a  listing  of  subroutine  BODES.  A  FORTRAN  subroutine  is  used  since  complex 
arithmetic  can’t  be  handled  in  ACSL.  A  dollar  sign  in  column  1  at  the  beginning  of  the  routine 
causes  the  ACSL  compiler  to  add  the  program  LABELED  COMMON  block  to  this 
subroutine.  Thus,  it  is  not  necessary  to  specify  arguments  in  a  subroutine  call. 

S,  G,  NUM  and  DEN  are  declared  complex  variables.  S  is  the  s  domain  variable  which  for 
steady  state  sine  waves  is 

S-jw 

or 

S=CMPLX  (0.0, W). 

G  is  the  complex  gain  of  the  transfer  function.  NUM  and  DEN  are  the  complex  values  of  the 
numerator  and  denominator,  respectively.  For  Z-transfer  functions 

Z-1  =  e”8^  =  e  ■  cos  wT  -  jsin  u>T 
or 

S=CMPLX  [cos  (W*DT),-SIN(W*DT)] 

G,  NUM  and  DEN  are  initialized  by 
G=CMPLX  (O-.p.) 


NUM=CMPLX  [A(l),0.] 


SUBROUTINE  BODES 

COMPLEX  S.  G.  NtjM*  OEN 
S>CMPLX(0.0*W) 

GO  TO  30 
ENTRY  BODEZ 

S*CMPLX(COS(W*DT>  *-SIN<W*Dt) ) 
30  G«CMPLX(0..0*> 

NUM*CMPLX (A ( 1 ) *0 • ) 

OEN*CMPLX(B(1)*0*> 

IF  (M«LE» 1 ) GO  TO  40 
00  10  1*2. M 
10  NUM»NUM«S*A(I) 

40  CONTINUE 

IF  (N*LEf«l>GO  TO  50 
00  20  1*2. N 
20  OEN*DEN*S.BCI) 

50  CONTINUE 
G*KG*NUM/DEN 
Y*AIMAG(G) 

X.PEAL(G) 

RETURN 

END 


Figure  4.  Listing  of  subroutine  bodes. 


DEN=CMPLX  [B(l),  0.] 


here  the  general  s  domain  transfer  function  looks  like 

A(l)  sm  +  A (2 ) sm_ 1  +  ...  +  A (m) s1  +  A (m+1) 

B (1)  sn  +  B<2)sn-1  +  ...  +  B(n)s1  +  B (n+1) 

and  the  general  Z  domain  transfer  function  appears  as 

A(l)Z_m  +  A(2)Z-m+1+  ...  +  A(m) Z_1  +  A (m+1) 

B  (1)  Z*"n  +  B  (2)  Z-n+1+  ...  +  A{n)  Z-1  +  A(n+1) 

The  latter  is  just  the  reverse  of  the  way  the  polynomials  of  (Z  * )  were  previously  written  in 
Section  2.  This  was  done  in  order  to  use  the  algorithm  originally  employed  for  theiS 
polynomials. 

Since  the  numerator  polynomial  can  be  written,  using  Horner’s  method,  as 

{„.  [(A(l)  s  +  A(2)).s+  A(3)J  s  +  A(4)  +  ...  °4) 

1  +  A  (n)  J  s  +  A  (n+1) 

then  a  DO  loop  can  be  used. 


DO  10  I  =  2,  M 


10  NUM  =  NUM  *S  +  A(I) 


where  M  is  defined  in  the  INITIAL  section  and  the  actual  value  used  is  set  in  the  run  time 
commands  (since  the  bilinear  transformation  and  impulse  invariant  methods  give  different 
order  numerator  polynomials).  Similarly  for  the  denominator  the  algorithm  is 


DO  20  I.«  2 ,N 


(15) 


20  DEN  =  DEN  *S  +  B(I) 
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When  the  numerator  or  denominator  consists  of  only  a  constant  term,  Af  —  1  or  N=  1 ,  then  its 
DO  loop  is  skipped. 

G=KG*NUM/DEN 

calculates  the  complex  transfer  function  gain  where  KG  is  an  optional  user  input  that  defaults 
to  1.0  if  not  specified.  The  real  and  imaginary  parts  of  G,  X  and  Y,  are  used  by  the  ACSL  main 
program.  Since  X  and  Y  are  named  in  the  ACSL  coded  portion  of  the  program  they  are  listed 
in  the  program  LABELED  COMMON  block. 

The  simulation  was  run  interactively  on  a  CYBER  74  using  ACSL  run  time  commands 
entered  on  a  Tektronix  console.  Hard  copy  plots  were  generated  along  with  line  printer 
listings.  Run  time  commands  that  were  used  every  session  were  put  on  a  temporary  mass 
storage  file  that  could  be  attached  as  a  user’s  local  file  to  the  Tektronix  terminal. 

Figure  5  is  a  list  consisting  of  interactive  session  setup  and  ACSL  run  time  commands.  After 
logging  onto  the  CYBER  74  INTERCOM  system  via  the  Tektronix  terminal  the  following  set 
up  statements  are  entered: 

CONNECT,  OUTPUT 

specifies  that  the  system  output  file,  default  name  OUTPUT,  will  be  displayed  on  the 
Tektronix  screen. 

ETL,  170 

extends  the  execution  time  per  statement  entered  to  170  octal  seconds. 

ATTACH,  INPUT,  TEK,  ID=DDXXXH  attaches  the  highest  cycle  (2)  of  TEK  to  the 
Tektronix.  INPUT  is  the  default  file  name  for  data  input  to  the  local  program  executing  in  the 
system. 

ATTACH,  LGOB,  TEK,  ID=DDXXXH,  CY—  1  attaches  the  simulation  absolute  binary 
file  which  was  previously  compiled  in  a  batch  job  via  fhe  system  control  cards  in  Figure  6. 

After  this  initial  preparation  LGOB  is  executed  by  entering 


LGOB. 


i 


CONNECT. OUTPUT 
ETL.170 

ATTACH. INPUT, TEK.ID*DOXXXH 
ATTACH » LGOB, TEK » ID*DDXXXH,CY*1 
L60B 


SET  TITLE <6>*"-7  HZ*  OTslONS" 

IT 

SET  OT*0 #020. TITLE (6) *M-7  HZ.  DT=20MS"SIT 


* 

S 


STOP 

BATCH. PPINT. PRINT. 3D.HDMD0 


Figure  5.  Listing  of  interactive  session  setup  and  ACSL  run  time  commands. 


HONDO, CN77000# 

ACCT • # • 

ATTACH. NACFIL.OCHACFIL ♦ ID-OCACSLSYS. 
ATTACH, ACSL. OCACSL.IOsOCACSLSTS# 

ACSL. 

FTN ( I*C0MPILE,B*2) 
request.lgob,*pf# 

NAP, OFF# 

ATTACH. ACSLLIB .OCACSLL IR  * I0*DCACSLSYS# 
ATTACH.GRAPH»TEKTRONlX40U,IDsWTPLOT.CY*2# 
LDSET .L IB*GRAPH, SUBST*ZZOtAW-TEKPLT  # 
LOSET.LlB»ifcCSLLIB.PPESET=INOEF. 

LOAO.LGO# 

NOGO.LGOB# 

RETURN. ACSL, NACFIL.ACSLLIS# 

CATALOG, LGOB.TEK, IO*OOXXXH. 

EXIT. 


Figure  6.  CYBER  74  control  cards  to  generate  an  absolute  binary  file,  LGOB,  of  the 
program  in  Figure  3. 


mm, 


■-  v om-coot ' 


The  ACSL  Executive  reads  the  ACSL  run  time  commands  from  logical  unit  number  5  (by 
default)  which  is  equivalent  to  reading  the  file  INPUT,  the  contents  of  which  are  shown  in 
Figure  7.  Referring  to  Figure  7, 

SET  PRN=9 

tells  the  simulation  that  the  line  printer  data  will  be  placed  on  file  PRINT  (logical  unit  number 
9)  instead  of  the  default  file  name  OUTPUT  (logical  unit  number  6).  After  termination  of  the 
simulation  the  PRINT  file  is  batched  to  a  line  printer. 

SET  PRNPLT=. FALSE.,  CALPLT=.TRUE.,  TTLCPL=.TRUE. 

replaces  the  default  printer  plot  routines  with  the  Tektronix  interface  routines  with  allowance 
made  for  plot  titles. 

SET  TITLE  (1)=  “2  JULY  IT 

puts  the  date  in  the  first  word  of  the  plot  title  array. 

PREPAR  T,  LOGW,  W,  WCPS,  GAIN,  PHASE,  MAG,  X,  Y  specifies  the  variables  to  be 
recorded  for  line  printer  listing. 

PROCED  ST 

SET  CNTNUS=.TRUE.,  FTSPLT=.TRUE.,  M=l,  N=4,  BILNR=.FALSE„ 
NVRNT=.FALSE. 

END 

defines  a  PROCED  which,  when  invoked,  picks  the  s  domain  version  of  the  Butterworth  filter 
for  Bode  plotting,  sets  the  order  of  the  numerator  and  denominator  polynomials  and  sets 
FTSPLT,  the  logical  for  suppressing  fly  back  tracing,  to  TRUE.  The  intent  here  is  to  overlay 
plot  either  Z-transform  method  with  the  s  domain  Bode  plot.  Two  runs,  one  for  the  s 
domain  and  one  for  the  Z,  are  to  be  made  without  rewinding  the  data  file.  So  when  a  plot  is 
made  the  overlay  is  also  automatically  generated.  The  FTSPLT  variable  set  TRUE  ensures 
that  a  line  connecting  the  end  point  of  the  first  plot  to  the  first  point  of  the  second  plot  will  not 
appear. 


SET 

>v 

ST 

•nV 

HO 

& 

ZB 

HI 

GO 

oooooooooooooooooooooo 

SET  PRN-9 

SET  PRNPLT*. FALSE.*  CALPLT».*RUE.*  TTLCPL*.TRUE. 

SET  TITLE <1)*"2  JUL  79“ 

PREPAR  T*  LOG*.  Vt  WCPS*  GAIN*  PHASE*  HAG*  X.  Y 
PROCED  ST 

SET  CNTNUS».TRUE.*  FTSPLT«.TRUE..  M«|*  N*4.  BILNR. •FALSE.*  NVRNT*. FALSE, 
END 

PROCED  ZN 

SET  CNTNUS* .FALSE • »  H»3*  N*A»  BILNR>JFALSE. *  NVRNT*.TRUE. 

ENO 

PROCED  ZB 

SET  CNTNUS*. FALSE.*  M»4.  ft«4*  BIlNR* JTRUE. •  NVRNT*.FALSE. 

ENO 

PROCED  GO 
START 

PRINT  "ALL" 

PLOT  “XAXIS"»LOGW*“XLO“*O.S*mXHIm«2.5 
DISPLY  MBCPS*  A*  B 
SET  NP*ITG*.FALSE. 

END 

PROCED  HO 
START 

DISPLY  WBCPS*  A*  B 
SET  NRNITG*.TRUE. 

END 

PROCED  IT 


PLOT  GAIN  S  SPARE 
PLOT  PHASE  S  SPARE 

SET  TITLE l2)««C0NT  VS  INV  DISCRETE  BUTTERVORTH  FILf" 

ST 

HO 

ZN 

GO 


PLOT  GAIN  f  SPARE 

W 

PL 01  PHASE  S  SPARE 

END  S  "IT" 

v 

SET  CMO-OIS 

Figure  7.  Tha  ACSL  standard  sat  of  run  tlmacommanda  residing  on  lodal  (No  INPUT. 


PROCED  ZN 


SET  CNTNUS=.F.,  M=3,  N=4,  BILNR=.F.,i  NVRNT=. T. 

END 

defines  a  PROCED  which  is  used  to  select  the  impulse  invariant  Z-transform  for  Bode 
plotting. 

PROCED  ZB  (16) 

SET  CNTNUS=.F.,M=4,  N=4,  BILNR=.T.,  NVRN=.F. 

END 

defines  a  PROCED  that  causes  the  bilinear  Z-transform  to  be  used  for  Bode  plotting. 
PROCED  GO 
START 
P/IINT  ‘ALL’ 

PLOT  ‘XAXIS  — LOGw-  ‘XLO’  =0.0,  ‘XHr=2.5 
DISPLY  WBCPS,  A,  B 
SET  NRWITG=. F. 

END 

defines  a  PROCED  that  executes  the  simulation,  prints  data  as  described  by  the  PREPAR 
statement,  describes  the  abscissa  for  the  Bode  plots  and  displays  the  break  frequency  and 
coefficients  for  the  transfer  function  numerator  and  denominator  polynomials.  Setting 
NRWITG  to  FALSE  causes  the  file  storing  the  PREPAR  list  variables  not  to  rewind  between 
runs  thereby  accumulating  all  of  the  data  from  the  runs  as  a  unit  to  be  printed  and  plotted 
together. 
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PROCED  HO 


START 


DISPLY  WBCPS,  A,  B 


SET  NRWITG=.T. 


is  a  PROCED  that  functions  much  the  same  as  ‘PROCED  GO’ but  it  turns  the  NRW1TG  flag 
to  TRUE. 

PROCED  IT 

SET  TITLE  (Z)=“CONT  VS  BILNR  DISCRETE  BUTTERWORTH  FILT” 

ST 

HO 

ZB 

GO 

PLOT  GAIN  $  SPARE 
PLOT  PHASE  $  SPARE 

SET  TITLE  (Z)=“CONT  VS  INV  DISCRETE  BUTTERWORTH  FILT” 

ST 

HO 

ZN 

GO 

PLOT  GAIN  $  SPARE 
PLOT  PHASE  $  SPARE 
END  $  “IT” 

is  a  PROCED  that,  when  invoked,  overlays  plots  for  the  bilinear  Z-transform  and  s  transfer 
function  and  for  the  impulse  invariant  Z-transform  and  s  transfer  function.  The  SPARE 
commands  following  the  PLOT  commands  automatically  generate  hard  copies  of  the  Bode 
plots. 

The  following  lists  the  code  for  FORTRAN  subroutine  SPARE, 

SUBROUTINE  SPARE 
CALL  HDCOPY 


I 


v.  *  ,  i.r, 


V.V.nV 


RETURN 

END 


HDCOPY  is  a  special  Tektronix  plot  software  routine  which  automatically  causes  the 
hardcopying  of  whatever  is  on  the  Tektronix  terminal  CRT.  The  name  SPARE  in  itself  is  a 
special  name  that  ACSL  recognizes.  SPARE  can  uniquely  be  used  for  the  execution  of  any 
routine  named  SPARE  during  run  time. 

SET  CMD=DIS 

indicates  to  the  ACSL  Executive  that  run  time  commands  will  now  be  input  from  the  DISPLY 
file  whose  logical  unit  number  is  DIS.  DIS  is  by  default  6,  which  is  also  the  number  for  the 
OUTPUT  file.  Thus  run  time  commands  are  now  expected  to  come  from  the  Tektronix 
terminal. 

Referring  to  Figure  5  the  following  commands  are  now  entered  on  the  Tektronix  display, 

SET  TITLE(6)=“-7  Hz,  DT=10MS" 

IT 

adds  to  the  plot  title  the  break  frequency  value  and  the  sampling  interval.  The  values  shown 
are  the  program  default  values.  PROCED  IT  is  then  invoked  and  automatically  generates  the 
plots  shown  in  Figurar8  through  //.  The  printer  listing  (not  shown)  shows  that  7  Hz  is 
approximately  the  -3  db  point  for  both  the  bilinear  Z-transform  and  the  s  transfer  function. 
However  the  drop  in  gain  afterwards  is  quite  dramatic  for  the  bilinear  method  as  shown  in 
Figure  8  where  GAIN  is  in  hundreds  of  db.  The  abscissa  range  in  frequency  value  from  1  to  50 
Hz  (LOGW=logio«u,  w  in  radians).  At  50  Hz  the  bilinear  method  has  a  gain  of -150db  and  the  s 
transfer  function  is  at  -51  db.  Figure  9  illustrates  the  good  phase  (in  degrees)  comparison 
between  the  two  methods.  The  maximum  deviation  is  approximately  18  degrees  at  50  Hz. 

As  shown  in  Figure  10  the  impulse  invariant  method  closely  follows  the  s  transfer  function. 
At  50  Hz  the  gain  of  the  impulse  invariant  transform  is  -56  db,  only  a  5  db  difference  from  the 
s-transfer  function.  However  in  Figure  11  one  can  see  the  phase  is  quite  a  bit  off  at  the  higher 
frequencies.  At  5Q  Hz  the  continuous  filter  has  a  phase  of  106  degrees  (actually  -254  degrees, 
the  ATANfc  routine’s  angle  limits  are  ±rr  radians)  and  the  digital  filter  is  at  1 75  degrees.  At  40 
Hz  there  is  approximately  a  10  degree  difference  in  phase. 

SET  DT=0.020,  TITLE  (6)=M-7  Hz,  DT=0OMS"$IT  sets  the  sampling  interval  at  20  ms 
and  generates  the  Bode  plots  in  Figures  12  through  15.  The  frequency  spectrum  repetition  is 
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Figures.  PhtNvmui  frequency  comparison  of  thas  domain  and  MttnearZ  domain 
filter  representation*  at  a  10  msec  sampling  interval. 
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Figure  1 1 .  Phase  vereue  frequency  comparison  of  the  s  domain  and  impulse 

invariant  Z  domain  filter  representations  at  a  10  msec  sampling! interval. 
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2  JUL  79  CONT  US  INU  DISCRETE  BUTTERUOR 
TH  FILT  -7  HZ-  DT-20MS 


Figure  15.  Phase  versus  frequency  comparison  of  the  s  domain  and  Impulse 
Invariant  Z  domain  filter  representations  at  a  20  msec  sampling  interval. 


just  starting  to  show  up  at  the  higher  frequencies.  In  Figure  /2the  bilinear  gai*  remains  below 
the  continuous  filter  out  to  31  Hz.  However,  in  Figure  13  the  phase  jumped  from  positive  to 
negative  at  25  Hz-the  reciprocal  of  twice  the  sampling  interval.  For  the  invariant  method  25 
Hz  is  the  cross  over  point  for  both  the  gain  and  phase  as  shown  in  Figures  14  and  15. 

DT  was  set  to  50  and  100  msec  with  the  results  shown  in  Figures  16  through  23.  Note  that  the 
bilinear  method  keeps  its  steep  roll  off  whereas  the  invariant  method  shows  the  frequency 
aliasing  phenomenon  that  destroys  its  roll  off  characteristics.  The  discontinuous  look  of  the 
bilinear  method  plots  is  caused  by  the  sparseness  of  data  points  at  the  higher  frequencies  due  to 
the  geometric  progression  method  of  generating  frequency  plot  points.  Figures  24  through  27 
are  plots  of  the  gain  and  phase  versus  frequency  (rad)  of  the  two  digital  filters  overlayed  onto 
the  continuous  filter  for  a  sampling  interval  of  100  msec.  Using  plain  frequency  for  the  abscissa 
brought  out  the  periodicity  of  the  digital  filters. 

STOP 

terminates  the  simulation.  The  PRINT  file,  now  attached  to  the  interactive  terminal  as  a  local 
file,  is  batched  for  line  printing  to  interactive  terminal  3D, 

BATCH,  PRINT,  PRINT,  3D,  HDMDD. 

The  sample  interval  value  DT  is  determined  by  the  frequency  content  of  the  input  signal 
since  input  signal  frequency  fold-over  can  occur  in  both  digital  filters  if  DT  is  made  too  large. 
However,  for  the  impulse  invariant  digital  filter  implementation,  the  actual  shape  of  the 
frequency  response  curve  itself  becomes  “aliased”  (as  mentioned  in  Section  2)  as 
1  cycle  approaches  the  break  frequency. 

2*DT  sec 

4.  TRANSIENT  RESPONSE  PROGRAM  AND  RESULTS  FOR  A 
THIRD  ORDER  BUTTERWORTH  FILTER  FOR  BOTH  THE  Z  AND 
s  DOMAIN  REPRESENTATIONS 

An  ACSL  program  was  written  to  compare  the  transient  response  of  the  two  digital  filter 
implementations  with  the  analytical  transient  response  of  the  continuous  version  of  the  filter. 
The  inputs  are  sine  waves  and  step  functions.  Also  programmed  was  a  standard  digital 
simulation  method  that  uses  the  RK2  numerical  integration  scheme  to  implement  the  M- 
method  canonical  form  of  the  continuous  filter  (as  discussed  in  Section  2).  The  numerical 
integration  scheme  will  be  compared  with  the  digital  filter  schemes.  It  would  also  be  desirable 
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Figure  19.  Phase  versus  frequency  comparison  of  the  s  domain  and  impulse 

invariant  Z  domain  filter  representations  at  a  50  msec  sampling  interval. 
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Flgur*  23.  Phase  versus  frequency  comparison  of  the  s  domain  and  impulse 

Invariant  Z  domain  filter  representations  at  a  1 00  msec  sampling  interval. 
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Figure  25.  Phase  versus  frequency  comparison  of  the  s  domain  and  bilinear  Z 
domain  filter  representations  at  a  100  msec  sampling  interval. 


INU  DISCRETE  BUTTERUOR 


to  compare  the  steadied-out  values  from  the  transient  response  with  those  predicted  by  the 
Bode  plot  program. 

As  discussed  in  Section  2  and  Appendix  A  the  digital  filters  were  implemented  in  macro 
ZXSFRM  which  uses  the  M-method  canonical  form.  The  digital  filters  are  also  implemented 
in  the  program  using  Bellman's  ‘Rectangular  Method'  canonical  form  for  comparison 
purposes.  Figure  28  shows  the  general  form  for  the  Rectangular  method. 

The  ACSL  transient  response  program  listing  is  in  Figure  29.  Macros  BTRWT1  and 
BTRWT2  are  the  Rectangular  method  implementation  of  the  impulse  invariant  and  the 
bilinear  methods  for  a  7  Hz  break  frequency  and  10  msec  sampling  interval,  respectively. 
Macros  CNVRNT  and  CB1LNR  were  discussed  in  Section  3.  They  generate  the  polynomial 
coefficient  values  for  the  two  digital  filters  as  implemented  in  the  general  purpose  Z-transform 
macro  ZXSFRM.  Parameters  are  defined  next  and  the  program  I/O  are  listed.  The  I/O  will 
be  defined  as  the  program  is  delineated.  In  the  prelNITIAL  section  global  constants  are 
defined.  These  are  parameters  such  as  computer  arithmetic  minimum  and  maximum,  RMN 
and  RMX;  jt;  the  conversion  factor  for  radians  to  degrees,  RADDEG;  logical  variable 
DUMP  which  controls  the  debug  printing  of  the  program  variables  for  the  final 
communication  interval;  and  the  stopping  time  for  a  simulation  run,  TSTP. 

S 

C1NTERVAL  CINT=0.050 
NSTEPS  NSTP=1 

sets  the  communication  interval  that  is,  the  interval  of  time  between  plot  points  or  line  printer 
outputs. 

In  the  INITIAL  section  the  parameters  for  this  particular  program  are  definer  TV'*  input 
amplitude  and  frequency  for  a  sine  wave  and  the  start  time  for  the  step  function  *  c  give 
values.  Macros  CNVRNT  and  CBILNR  define  the  numerator  and  denominator  coefficient 
U  and  V  for  the  impulse  invariant  method  and  R  and  S  for  the  bilinear  method.  Constant  gains 
for  the  analytical  response  of  the  continuous  filter  to  a  unity  gain  sine  wave,  Ka,  K./J,  K-y,  K6 
and  Kt,  are  calculated.  These  constants  correspond  to  A,  B,  C,  -D  2/  \/3  and  E/  wa  of  Equation 
(4)  in  Section  2. 

In  the  Dynamic  section  the  analytical  continuous  filter  response  is  calculated  according  to 
Equations  (3)  and  (4).  The  output  is  named  YANALF.  The  TERMT  macro  is  used  at  the  end 
of  the  DYNAMIC  section  to  stop  the  simulation  when  global  time  T  exceeds  the  value  in 
variable  TSTP. 


Figure  28.  Bellman’s  rectangular  method  canonical  form. 


PROGRAM  TRANSIENT  RESPONSE  FOR  CONTINUOUS  AND  DISCRETE  BUTTERWORTh 


MACRO  BTRMTl(BtA) 

MACRO  REDEFINE  X , Y*XIC*YIC» A1 * A2.A3.B1 *B2* I 
MACRO  RELABEL  SN©01*SN002 
ARRAY  XIC<3)  »YIC<3>  .X(3)  .¥(3) 

CONSTANT  YIC=3*0.0*XIC=3*0.0 

CONSTANT  Bl=3*1381538«E-2*B2=2.34235779E-2.Als*2.l34;>8885.A2=-l. 60402376 
CONSTANT  A 3= ♦0*414929794 
PROCEDURAL <B=A> 

IF (ZZFST (I) .LT.0»5)G0  TO  SN002 

DO  SN001  1=1*3 

CALL  ZZICS(YICd)  »  Y  ( I ) ) 

SNOO 1 • .CALL  ZZICS<X1C(I)«X<!>> 

B=Y(1)*A1*Y(2)*A2*Y<3)*A3*X(1)*B1*X(2>*B2 

Y(3)*Y(2> 

Y (2) *Y  (1 ) 

Y ( 1 ) *B 
X (2) *X  ( 1 ) 

X  ( 1 )  =*A 

SN002.. CONTINUE 
ENO 

MACRO  END 
MACRO  BTRMT2(B»A) 

MACRO  REDEFINE  X, Y*XIC*YIC.I *A1 *A2»A3,B1 .B2.B3.B4 
MACRO  RELABEL  SN001«SN002 
ARRAY  XIC(3) * YIC  <3 ) »X ( 3) » Y  (3 ) 

CONSTANT  YIC»3*0.0*XIC*3*#.0 
CONSTANT  Al=**2. 12664  ... 

*A2=-1. 59582*A3«0. 41 184*81=7. 167e-3*B?=2.1503E~2*b3=?.1503E-2.B4*7. 16... 
7E-3 

procedural <e*A) 

IF (ZZFST (I). LT. 0.5)60  TO  SN002 

DO  SN001  1=1*3 

CALL  ZZICS(YIC(I).Y(I)) 

SN001..CALL  ZZICS(XIC(I) *X(I>) 

B=Y(l)*Al.Y(2)*A2.Y(3)*A3*BlPA.B?«X(l).B3*X(2).B4*Xn) 

Y (3) =Y (2) 

Y (2) =Y  < 1 ) 

Y ( 1 ) =B 
X (3)*X (2) 

X (2) =X  < 1 ) 

X(1)=A 

SN002.. CONTINUE 
END 

MACRO  ENO 


Figure  29.  Listing  of  the  ACSL  transient  response  program. 
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MACRO  CNVRNT<U»V*WB*DT*KG> 
macro  redefine  k*e*c*s 

ARRAY  U(3).V(4> 

PROCEDURAL <U*V*WB*DT*KG) 

E=EXP(-WB*DT) 

C=SQRT  <E)  *COS  (SORT  (3* )  *0.5**B*DT) 

S=SORT  (E )  *S IN (SORT  ( 3  • )  *0#5*k(fl*DT )  /SORT  (3* ) 

U(1)=0.0 

U(2)=WB*(-C>S4E) 

IM  3) »WB*E*  < 1 *-C-S) 

V(1)31»0 
V  <2) =-2.*C-E 
V(3>**E*<1.*2**C> 

V (4 ) =-E**2 

K=(U<2>*U(3>)/<V<1)*V<2MV(3MV<4>> 

U(1)=U(1)*KG/K 

U(2)=U(2)*KG/K 

U (3) =U(3) *KG/K 

END 

MACRO  END 

MACRO  CBILNR(R*StWB*DT*KG) 

MACRO  REDEFINE  AL.DEN 
ARRAY  R  (4 ) tS  (4 ) 

PROCEDURAL (R* S*W8f OT*KG) 

AL=TAN(0.5«WB*DT) 
DEN=1.*2.*AL*2#*AL**24AL**3 
R(1)=KG*AL**3/DEN 
R (2) =KG*3.*AL**3/DEN 
R  (3) =R (?) 

R (4 ) =R ( 1 ) 

S<1)*1.0 

S(2)=(-3,-2.*AL^2.*AL«*?*3.«AL**3)/DEN 

S<3)=<3.-2.*AL-2.*AL**2*3.*AL**3)/DEN 

S(4)=(-1.42.*AL-2,*AL«*24»L**3)/0EN 

END 

MACRO  ENO 


Figure  29.  (Continued) 


MACRO  MACRO  ZXSFRM (P»Q) 

MACRO  REDEFINE  V.I»YIC*J#VP 

MACRO  RELABEL  SNOOl .SNOO?*SN003tSN004.SN005 

MACRO  ASSIGN  N 

MACRO  MULTIPLY  O 

MACRO  INCREMENT  Q<4) 

MACRO  01 VIDE  0(3) 

MACRO  IF<N*0>999 
MACRO  MULTIPLY  0 
MACRO  INCREMENT  Q<4) 

ARRAY  Y(N)«YIC(N> 

CONSTANT  YIC*N#0.0 
PROCEDURAL (P <1 > >P (2) «P ( 31 *P  (4> ) 

IF (ZZFSTtI) .LT. 0.5)60  TO  SN005 
DO  SN041  1-l.N 

SN001..CALL  ZZICS(YIC(I),Y(I>) 

YP-P  <  ?  > 

MACRO  DECREMENT  1 
DO  SN002  1*1 «N 
MACRO  INCREMENT  1 
SNo02..YP-YP-P<4> <I«1)*V<N-I» 

Y(N)*YP/P(4M1) 

VP-0.0 

MACRO  INCREMENT  1 
DO  SN003  1*1 «C  <3) 

SN003. • YP*VP«P  (3)  m«Y<N-I» 

P  < 1 1 -YP 

MACRO  DECREMENT  2 
DO  SN004  1*1 *N 
SN004. #Y ( I ) *Y ( I ♦ 1 ) 

SN005.. CONTINUE 
END 

MACRO  EXIT 

MACRO  999.. PRINT  NUMERATOR  GREATER  THAN  DENOMINATOR 
MACRO  ENO 

"PARAMETERS" 

"A-AMPLITUOE  OF  INPUT  SIGNAL  TO  FILTER" 

"MMTHO-LOGICAL*  IF  .T.  INOICATES  M-METHOO  IMPLEMENTATION  IS  TO  BE  USED" 
"STPSIN-LOGICAL  INDICATING  WHETMeR  INPUT  SINE  WAVE  C.F.)  OR  STEP  < . T . >  •• 
"TZ-TIME  AT  WHICH  STEP  TURNS  ON  -  SEC" 

"W-INPUT  SINE  WAVE  FREQUENCY  -  RAO/SEC" 

"WBCPS-F ILTER  BREAK  FREQUENCY  -  HZ" 

"WB-FILTER  BREAK  FREQUENCY  -  RAO/SEC" 

"WCPS-W  IN  HZ" 

"INPUTS" 

"MMTHD ♦  STPSINt  WBCPS.  A.  WCPS.  TSTP*  CINT*  MAXT 1 «  MAXT?*  MAXT3,  TZ" 
"OUTPUTS" 

"YINTEG.  U«  V*  R*  S,  YINV,  YBILNR.  YANALS*  YANALF*  YANALRt  XZFERt  ... 
XINTEG" 


Figure  29.  (Continued) 


h - GLOBAL  CONSTANTS 

LOGICAL  DUMP 


CONSTANT  RMN  «  1.0E-30*  RMX  »  1.0E30  •  RAD0EG*57 • 2957795 
.  PI»3. 14159265*  l»UMP*JF ALSE.  *  TSTP  =  3.0 
TWOPI  »  2.0«PI 
CINTERVAL  CINT  >  0.650 
NSTEPS  NSTP  «  1 


INITIAL 

- - - - SET  CONSTANTS  FOR  INPUT  SIGNAL  AND 

DESIRED  BREAK  PT  FREOUENCY  OF  FILTERS  •• 

LOGICAL  STPS1N 

CONSTANT  STPSIN  *  .TRUE. 

CONSTANT  WCPS  *  50,  .  A  *  1)0  ,  TZ  =  0.049999999 

W  *  WCPS*TWOPI 

CONSTANT  WBCPS  >7.0 
WB  *  WBCPS«TWOPI 

" - - - INVARIANT  FILTER  COFFFS  FOR  3RD  ORDER  ... 

BUTTERWORTH  •• 

CNVRNT4U*  V»WB*  MAXT2,  1.0) 

" - - - BILINEAR  FILTER  COEFFS  FOR  3RD  ORDER  ... 


butterworth  " 

CB1LNR(R.  S*WB*  MAXT2*,  1.0) 

- -----constants  for  analytical  response  of  ... 

FILTER  to  a  sine  wave  w/unity  gain  « 

KAL  *  WB**3«(W**3  -  2.0#WB**2*W)/(W«*6  ♦  WB«*6) 

KBE  *  WB**4« (-2«0*W#*2  ♦  WB«*2)/(W»*6  ♦  W9*«6) 


KG A  *  W8*W/(WB»*2  ♦  W**2> 

KDL  »-1.1547*WB*(W**3  -  WB**2*W) / (W«*4  -  W**2*WB«*2  ♦  WB**4> 
KEP  *  WB«*3*W/ (•••A  -  W«*2*WB**2  ♦  *8**4! 


ENO  S  "INITIAL" 
DYNAMIC 


yanalf  * 


- calculate  the  analytical  value  for  the  ... 

BUTTERWORTH  FILTER  FOR  THE  FAST  DERIVATIVE  SECT" 
RSW (STPSIN*  RSW (T.GE6TZ •  STEP(TZ)  -  EXP (-WB* (T-TZ) )  ... 

-1.1547*EXP(-0.5*WB*(T-TZ)>  ... 

•SIN  (0.86*025*WB«  CT-TZ))  .... 

•  KAL*COS(W*T)  ♦*KBE#SINIW*T)  ♦  KGA*EXP  III 
<-WB*T)  -  KDL*EXP(-0.5*WB«T)*SIN  ... 

(0.B66025«WB«T  -  1.947198)  *  KEP*EXP  ... 
<-0«5*WB*T)*SIN(WB*T)  ... 


DERIVATIVE  FAST 

MAXTERVAL  MAXT1  *  O.Old 
ALGORITHM  IALG  «  4 


TF  »  1NTEGO.O*  6.0) 

" - ... - - - - —CALCULATE  INPUT  TO  XSFER  FN 

XINTEG  «  RSW (STPSIN*  A*STEP(TZ)*  A«SIN(W*TF»> 

- - - ENUMERATE  XSFER  FN  COEFFICIENTS 

REAL  P<1)  *  0(4) 


CONSTANT  P  ■  6.508I2E4*  Q  ■  1«0  *  87.9646  • 

*  3868.89  *  8.50812E4 

» - - - - - NUMERICALLY  INTEGRATE  FOR  FILTER  OUlPUT" 

TRAN ( YINTEG*0»  3*  P*  0*  XlNTEO) 

ENO  S  "DERIVATIVE" 


Figure  29.  (Continued) 


DERIVATIVE  slo* 

- - - - - MA*T2  REPRESENTS  SAMPLE  INTERVAL  LENGTH 

MAXTERVAL  HAXT2  *  0.018 
ALGORITHM  IALG  =  4 


TS  *  INTEG(1.0,  0*0) 

- - CALCULATE  FILTER  OUTPUT  USING  THE  INV  . 

and  BILINEAR  Z-TRANSFORM  METHODS" 

- - - -CHOOSE  M-METHOD  FOR  IMPLEMENTING  FILTERS 

LOGICAL  MMTHD 

CONSTANT  HMTHOs.TRUE. 

« - .——INVARIANT  filter 

YINV  *  RSW(MMTHD.  ZXSERM (XZFERt  Ut  V) »  BTBWT1 (XZFER) ) 

n - - - - - bilinear  filter 

YSILNR  «  RSW (MMTHD*  ZXSFRm (XZFER*  R»  S) ♦  BTRWT2 (XZFER) ) 


PROCEDURAL  <  Y ANALS . XZFER* ) 

IF (ZZFST (TZ) .LT.0.5)G0  TO  SN082 

m - - - --CALCULATE  INPUT  TO  Z  TRANSFORM  FILTERS 

XZFER  *  RSW (STPSIN*  A*STeP(TZ>*  A«SIN(W*TS>) 

- - CALCULATE  THE  ANALYTtCAL  VALUE  FOR  THE  • 

butterworth  filter 

YANALS  *  RSW (STPSIN*  RSW(TS.GE.TZ*  STEP(TZ)  -  EXP(-WB*(TS-TZ)>. 

-1.1547*E*P(-0.5«WB*(TS-TZ>). 
•S|N<0.8«6625*WB*(TS-TZ>)  • 

*  (<i)  • 

•  KAL*COS <W*TS)  ♦  KBE«SIM(W«TS)  ♦  K6A-EXP  * 
«-WB*TS)  -  KDL*EXP(-0.6*W8*TS)*SIN  • 

(0.8660Z5*WB*TS  -  1.049198)  ♦  KEP*EXP  . 

<-o;5»wb«ts)»sin(vb*ts)  • 

) 

SN002.. CONTINUE 

ENDS"PROCEDm 


END  S  "DERIVATIVE" 
DERIVATIVE  RECORO 


- - -----SAMPLE  FAST  ENOUGH  TO  SEE  THE  STAIRCASE. 

IN  YINV  AND  YBILNRm 
MAXTERVAL  MAXT3*0.0920000I 
- - — - — UNSYNCHRONIZED  EULER" 


ii 


•• 


ALGORITHM  IALG*3 

- --establish  local  time  FOR  THIS 

DERIVATIVE  SECTION" 

TR  *  INTEGd.Ot  0#«> 

- -----CALCULATE  THE  ANALYTICAL  VALUE  )0R  THE 

butterworth  filter 

YANALR  ■  RSW (STPSIN*  RSW (TR.GE.TZ,  STEP(TZ)  -  fXP(-WB*<TR-TZ) ) 

-1 . 1547*EXP (-0.5«WB* (TR-TZ) ) 


* 


*S|N(0.8060a5»WB«(TR-TZ)) 

_  *  0.0) 

»  KAL*COS(W«TR)  ♦  KBE*SIR(W«TR)  ♦  K6A«EXP 
(-WB.TR)  -  NOL*EXP f-0.5*W8*TR) *SIN 
<0.86602S*WB*TR  -  1.049198)  *  EEPwEXP 

<-o;5#wb«tr)«sin(wb*tr> 


) 


- - — ...... — — — —RECORO  THE  DATA" 

LOGICAL  RECORO 


CONSTANT  RECORO*. FM.SE . 
IF (RECORO) CALL  ZZLOO 
ENO  S  "DERIVATIVE" 


TERNT(TfcGT.TSTP) 
ENO  S  "DYNAMIC* 


Figure  29.  (Continued) 


TERMINAL 

IF (DUMP) CALL  DEBUG 
END  $  "TERMINAL" 

END  $  "PROGRAM" 


1  Figure  29.  (Concluded) 


There  are  three  DERIVATIVE  sections:  DERIVATIVES  FAST,  SLOW  and  RECORD. 
DERIVATIVE  FAST  contains  the  numerical  integration  of  the  continuous  filter  as 
implemented  in  the  ACSL  TRAN  macro.  MAXT1  is  the  maximum  integration  step  size  as 
declared  by  MAXTERVAL;  the  RK2  algorithm  is  selected  through  the  ALGORITHM 
declarative  statement;  and  TF,  the  local  time  for  this  DERIVATIVE  section,  is  generated  by 
integrating  a  unit  step  function.  XINTEG,  the  input  to  TRAN,  is  either  a  step  or  sine  wave 
depending  on  whether  STPSIN  is  TRUE  or  FALSE,  respectively.  YINTEG,  the  output  of 
TRAN,  is  the  Butterworth  filter  output.  P  and  Q  are  the  numerator  and  denominator 
coefficients  of  the  continuous  filter  as  configured  in  Equation  (5)  of  Section  2.  TRAN  is  a 
standard  ACSL  system  macro  [1]. 

In  DERIVATIVE  SLOW  the  digital  filter  responses  are  calculated.  Instead  of  specifying 
[integration  step  size,  MAXT2  has  a  value  that  is  equivalent  to  the  sampling  interval  for  the 
digital  filters.  Setting  IALG  to  4  specifies  the  RK2  integration  algorithm  is  to  be  used  in  this 
DERIVATIVE  section.  The  RK2  algorithm  has  no  use  in  this  section  other  than  to  verify  that 
ZXSFRM  is  working  properly  in  a  DERIVATIVE  section  employing  a  multi-step  integration 
algorithm.  TS  is  the  local  time  in  this  DERIVATIVE  section.  The  logical  variable  MMTHD  is 
used  to  select  one  of  two  canonical  forms,  the  M-method  or  the  Rectangular  method.  YINV  is 
the  digital  filter  output  for  the  impulse  invariant  method  and  YBILNR  is  the  digital  filter 
output  for  the  bilinear  method.  XZFER  is  the  digital  filter  input  signal.  YANALS,  the 
analytical  continuous  filter  response  is  calculated  according  to  Equations  (3)  and  (4)  of 
,  Section  2. 

IF  (ZZFST(TZ). LT.0. 5)  Go  to  SN002 

9 

SN002... CONTINUE 

jumps  over  the  YANALS  and  XZFER  calculations  unless  the  ACSL  Executive  is  executing 
the  first  minor  step  of  the  multistep  integration  algorithm.  This  code  is  implemented  to  ensure 
*that  YANALS  and  YINV  and  YBILNR  and  XZFER  are  synchronized  in  time,  in  order  to 
compare  the  phases  of  the  transient  response  data.  It  is  to  be  noted  that  the  positions  in  time  of 
YINV  and  YBILNR  correspond  to  a  zero  delay  in  the  time  between  when  the  input  is  received 
in  the  filters  and  when  the  output  is  calculated  and  sent  out.  This  causes  no  problems  inside  a 
digital  computer  but  if  these  values  were  D  to  A’d  to  a  continuous  system  thei'e  is  some 
inherent  lag  which  probably  should  be  represented.  Since  the  digital  filter  input  and  output 
equations  are  calculated  only  at  the  front  of  a  major  integration  step  they  will  appear  in  line  - 
.printer  listings  and  plots  to  have  a  one  time  step  delay.  However,  values  for  YINTEG, 
YANALF  and  XINTEG  do  match  correctly  with  output  time  values. 


DERIVATIVE  RECORD  is  used  specifically  for  recording  data  at  a  faster  rate  (without 
changing  CINT)  in  order  to  show  the  staircase  nature  of  the  digital  filters  as  seen  by  a 
continuous  system.  As  mentioned  above  no  time  delays  due  to  the  time  used  for  the  actual 
calculation  of  the  filters,  algorithms  are  modeled.  Time  skewing  between  digital  filter  inputs  ' 
and  outputs  is  due  solely  to  the  intrinsic  features  of  the  filter  algorithms  themselves.  Y  AN  ALR 
is  the  analytical  filter  output  calculated  in  DERIVATIVE  RECORD  which  is  used  in  plotting 
as  a  comparison  with  the  staircase  digital  filter  outputs.  MAXT3  is  approximately  2  msec 
which  will  be  the  data  recording  interval.  It  is  made  slightly  greater  than  2  msec  to  make  sure 
the  digital  filters  get  iterated  on  the  first  minor  step  in  DERIVATIVE  SLOW  before  a  data 
logging  action  takes  place.  IALG  is  made  3  which  corresponds  to  the  synchronized  Euler 
integration  (one  DERIVATIVE  calculation  per  integration  step)  for  this  section.  If  RECORD 
is  TRUE,  a  data  logging  action  takes  place  through  the  call  to  ZZLOG 

The  simulation  was  run  interactively  on  a  CYBER  74  using  ACSL  run  time  commands 
entered  on  a  Tektronix  console.  Hard  copy  plots  were  generated  along  line  printer  listings. 

The  interactive  session  setup  commands  are  given  in  Section  3.  ACSL  run  time  commands 
that  were  used  every  session  were  put  on  a  file.  These  commands  are  listed  in  Figure  30. 
Discussion  of  this  run  time  command  file  will  center  on  the  PROCED’s  as  the  other 
commands  are  the  same  as  those  discussed  in  Section  3. 

ACTION  ‘VARMJ.039,  ‘VAL^O,  ‘LOC’=NDBUG 

causes  a  debug  printing  starting  at  0.039  secs  or  thereafter  with  NDBUG=20.  ‘PROCED  GO’ 
is  much  the  same  as  that  used  in  Section  3  with  the  addition  of  an  ACTION  ‘CLEAR’ 
statement  since  only  the  first  run  is  to  have  a  debug  printing. 

‘PROCED  MSTEP’  and  ‘PROCED  RSTEP’  configure  the  simulation  for  running  the  M- 
method  and  Rectangular  method  implementations  of  the  digital  filters,  respectively,  for  a  step 
.input.  The  analytical  and  numerical  integration  values  of  the  continuous  filter  are  always 
calculated.  Similarly  MS1N05  through  MSINSO  run  the  simulation  with  S  Hz  through  SO  Hz 
sine  waves  for  the  M-method  implementation.  RSIN40  sets  the  simulation  up  with  a  40  Hz 
sine  wave  input  and  the  rectangular  implementation  method  for  the  digital  filters. 

7  Hz  is  the  only  filter  break  frequency  used  in  the  following  runs.  Figure  31  through  33  are 
generated  by  the  following  run  time  commands. 


MSTEP 

GO 


Set  PRN-9 

SET  PRNPLTa.FALSE.*  CALPLT*.TRUE. .  TTLCPL=.TRUE. 

PREP AR  T , Y ANALS . YANALF » Y INTEG » Y INV ♦ YB I LNR ♦ XZFER • X INTFG ♦ TS . TF 
PREPAR  T.TR.YANALR 
SET  DUMP-.T. 

ACTION  "VAR"a0.039*"VAL"=20*"L0C"»NDBUG 
SET  TITLE C1)“M15  AUG  79-3RD  ORD  BTRNTH-7HZ" 

PROCEO  GO 
START 

ACTION  "CLEAR” 

PRINT  "ALL" 

DISPLY  U.V.R.S  S  "DISPLAY  DIG  FIlT  COEFFS" 

PLOT  "XHI"*T 
END 

PROCED  MSTEP 

SET  TITLE (A) »"STEP  INPUT  M-METHoO  " 

SET  MMTHDa.TRUE. 

SET  STPSIN-.TRUE.*  TSTP=0.5*  CINT=0.010 
END 

PROCEO  RSTEP 

SET  TITLE (A) »"STEP  INPUT  RECTANGULAR  METHOD  " 

SET  MMTHD». FALSE. 

SET  STPSIN-.TRUE**  TSTP*0.5.  CINTs0.010 
ENO 

PROCEO  MS I NOS 

SET  TITLE <4>-"SINE  INPUT  BHZ.  M-mETHOD  " 

SET  MMTHD-.TRUE. 

SET  STPS1N-.FALSE.*  TSTP-t.O*  CInT-OIOIO*  WCPS-5.0 
END 

PROCEO  MSINIO 

SET  TITLE <4)««SINE  INPUT  IOHZ*  M-METHOD  « 

SET  MMTHD».TRUE. 

SET  STPS1N-.FALSE.*  TSTP*0.5.  CInT«0*010,  WCPSMO.O 
END 

PROCEO  MSIN20 

SET  TITLE <4> -"SINE  INPUT  20HZ*  M-METHOD  » 

SET  MMTHO-.TRUE. 

SET  STPSIN*. FALSE.*  TSTP«0.4*  ClNT-OfcOlO*  WCPS»20.0 
-NO 

PROCEO  MSIN40 

SET  TITLE <4) -"SINE  INPUT  *OHZ*  M-METHOD  « 

SET  MMTHO-.TRUE. 

SET  STPSINa.FAi.SE.*  TSTP-O.A*  CInT*0«010*  WCPS-40.0 
ENO 

PROCED  RSINAO 

SET  TITLE <4>»"$INE  INPUT  4#HZ*  RECTANGULAR  METHOO  " 

SET  MMTHDv.FALSE. 

SET  STPSIN-. FALSE..  TSTP-0.4.  ClNT-0«010*  MCPS-4O.0 
END 

PROCEO  MSIN47 

SET  TITLE (4>«"SINE  INPUT  47MZ*  M-METHOD  " 

SET  MMTHO«.T. 

SET  STPSIN-. F*.  TSTP-0.5*  CINT-O.OIO*  NCPS-47. 

ENO 

PROCED  MSIN50 

SET  TITLE i4) -"SINE  INPUT  SOHZ*  M-METHOD  « 

SET  MMTHO-.TRUE. 

SET  STPSIN-.FALSE.*  TSTP-0.4.  CInT«0«010*  NCPS-50.0 
8?  CMD-DIS 


Figure  30.  A  standard  set  of  ACSL  run  time  commands  that  were  stored  on  file 
INPUT. 


Figure  32.  Plot  of  filter  Input  and  resulting  bilinear  digital  and  analytical  filter 
reap  once  using  a  10  msec  sampling  interval. 


PLOT  YBILNR.  YINV 

Pl.O'l  YANALS,  YBILNR,  XZFER,  ‘HL=2. 

PLOT  YANALF,  YINTEG,  XINTEG,  ‘HI’=2. 

Figure  31  shows  that  the  two  z-transform  methods  give  essentially  the  same  transient  response 
to  a  unit  step.  Figure  32  is  a  comparison  plot  between  the  analytical  and  bilinear  filter  values. 
The  bilinear  filter  leads  the  analytical  Filter  response  by  about  5  msec.  Figure  33  shows  how 
well  the  analytical  and  numerical  integration  responses  compare.  YINTEG  is  a  bit  low  at  its 
peak  amplitude  but  it  has  excellent  phase  matching.  Special  effort  was  made  to  synchronize  all 
the  filter  implementations  with  their  analytical  responses.  In  the  program,  TZ  was  made 
slightly  less  than  0.05  sec  to  ensure  that  the  STEP  macro  went  high  at  0.05  sec.  So  observed 
phase  differences  should  be  accurate  and  not  a  function  of  program  calculations  and  1/0 
functions. 

RSTEP  $  GO 

PLOT  YBILNR,  YINV 

PLOT  YANALS,  YBILNR,  XZFER,  ‘HI ’=2 

generates  data  using  the  Rectangular  method  canonical  form.  Plots  are  shown  in  Figures  34 
and  35.  The  results  are  identical  to  the  above  where  the  M-method  canonical  form  was  used  to 
implement  the  digital  filters.  Thus  both  canonical  forms  give,  as  they  should,  the  same  digital 
filter  response. 

SET  RECORD=.T.  $  GO  $  PLOT  YANALR.  YBILNR 

reruns  the  simulation  and  causes  the  data  recording  interval  to  drop  to  2  msec  which  is  used  to 
bring  out  the  staircase  nature  of  Y BI LN R  as  would  be  seen  by  a  continuous  system.  YANALR 
is  simply  the  analytical  filter  value  updated  every  2  msec.  Although  the  digital  filter  actually 
leads  the  analytical  filter  by  about  one-half  of  the  digital  filter  sampling  interval.  Figure  3$  a 
continuous  system  would  see  a  staircase  function  that  is  practically  in  phase  with  the  analytical 
filter  response.  Of  course,  the  delay  in  YBILNR,  due  to  the  finite  time  needed  to  calculate  the 
filter  algorithm,  has  been  neglected.  To  a  continuous  system  a  discrete  filter  signal  in  phase 
with  the  continuous  analytical  filter  signal  would  appear  lagging  in  phase  by  one-half  of  a 
sampling  interval  as  shown  in  Figure  37. 


MSIN05  $  SET  RECORD  =.F.$GO 
PLOT  YBILNR,  YINV 


VV  J  J-  «* 


{  r  I 


17  AUG  79 -3RD  ORD  BTRUTH-7HZ  STEP  INPUT 
RECTANGULAR  METHOD 


Figure  35.  Plot  of  filter  input  and  corresponding  responses  from  the  bilinear  digital 
filter  using  the  rectangular  canonical  form  and  the  analytical  filter  with  a 
10  msec  sampling  interval. 


WITT 


PLOT  YANALS,  YBILNR,  XZFER 
PLOT  YANALF,  YINTEG,  XINTEG 


runs  the  simulation  with  the  5  Hz  sine  wave  input  and  with  the  M -method  canonical  form  for 
the  digital  Filters.  Figures  38  through  40  are  generated  from  these  commands.  Figure 38  shows 
that  the  bilinear  and  impulse  invariant  digital  filters  compare  quite  well.  Figure 39  shows  that, 
except  for  an  initial  lead  transient  of  about  5  or  6  msec  in  the  bilinear  filter,  the  analytical  and 
bilinear  responses  are  similar  in  amplitude  and  phase.  Table  1  lists  the  amplitudes  and  phases 
for  all  four  filter  representations.  The  analytical  response  is  calculated  in  two  ways: 

•  From  the  analytical  expression  as  t— 00  and 

•  From  calculations  based  on  the  transient  response  data.  The  digital  integration 
method’s  amplitude  and  phase  are  also  calculated  based  on  transient  response  data.  Bode  plot 
data  for  phase  and  amplitude  are  listed  in  the  table.  All  phase  and  amplitude  calculations  are 
based  on  expressions  developed  in  Appendix  B. 

All  three  sources  of  data  for  the  analytical  filter’s  amplitude  and  phase  give  essentially  the 
same  results.  The  closeness  of  the  two  digital  filters’ amplitudes  and  phases,  as  calculated  from 
transient  response  data,  to  both  their  Bode  plot  data  and  the  analytical  filter  response  can  be 
seen  in  the  table. 

Similarly  the  numerical  integration  form  of  filter  implementation  has  an  output  whose 
amplitude  and  phase  match  those  of  the  analytical  filter  quite  well  as  shown  in  Figure  40  and 
Table  1.  RECORD  is  set  to  TRUE  in  order  to  generate  Figure  41  where  the  discrete  bilinear 
output  is  compared  to  the  analytical  filter  as  calculated  every  2  msec. 

Next  a  10  Hz  sine  wave  was  input  to  the  filters.  Figure  42  indicates  that  the  two  digital  filters 
have  similar  phases  but  the  amplitudes  are  now  starting  to  differ  because  of  the  faster  roll-off 
of  the  bilinear  filter.  Figure  43  compares  the  invariant  digital  filter  with  the  analytical 
response.  The  closeness  with  which  the  digital  filters  match  their  Bode  amplitude  and  phase  as 
well  as  those  of  the  analytical  response,  can  be  seen  in  Table  1. 

Figure  44  compares  the  numerical  integration  scheme  with  the  continuous  analytical  filter 
The  phase  match  is  quite  good  but  the  amplitudes  do  not  correspond  as  well.  Table  1  shows  the 
phase  of  the  numerical  integration  scheme  has  about  4  degrees  of  lead  over  the  analytical  filter 
and  has  approximately  16  percent  more  amplitude. 
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Figure  40.  Plot  of  filter  Input  and  filter  outputs  for  the  Integration  method  and  the 
analytical  expression  with  a  10  msec  Integration  step  size. 


73 


(v)  sibNyx 


Figure  43.  Plot  of  filter  Input  and  filter  outputs  for  the  Impulse  invariant  method  and 
analytical  expression  using  a  10  msec  sampling  Interval. 
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Figure  44.  Plot  of  filter  input  and  filter  outputs  for  the  numerical  method  and 
analytical  expression  with  a  10  msec  integration  step  size. 
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A  20  Hz  sine  wave  was  used  to  drive  the  filters.  Figure  45  is  a  comparison  plot  of  the  two 
digital  filters.  Again  their  phases  match  and  their  amplitudes  compare  with  their  respective 
Bode  data  as  shown  in  Table  l.  In  Figure  46  the  invariant  digital  filter  is  compared  with  the 
analytical  response.  They  are  essentially  the  same  after  the  first  transient  oscillation.  Figure  47 
and  Table  1  show  that  although  the  integration  method  has  58  percent  more  amplitude  than 
the  analytical  solution  it  lags  the  analytical  phase  by  only  7  degrees. 

A  40  Hz  sine  wave  input  resulted  in  digital  filter  outputs  as  shown  in  Figure  48.  Figure  49 
compares  the  impulse  invariant  and  analytical  responses.  Numerical  values  for  amplitudes  and 
phases  are  shown  in  Table  1.  Notice  how  the  sampled  input  sine  wave  bears  little  resemblance 
to  its  continuous  counterpart.  The  sampled  input  essentially  has  two  peaks  that  alternate  in 
time  and  the  filter  outputs  have  the  same  look.  The  filters’  amplitudes  and  phases  were 
calculated  from  the  larger  of  the  two  peaks  (for  example,  points  (3)  and  (4)  in  Figure  48  since 
using  low  peak  data  (for  example,  points  (1)  and  (2)  in  Figure  48)  gave  lower  amplitudes. 

Using  points  (3)  and  (4),  Equations  (B-6)  and  (B-7)  in  Appendix  B  give  the  amplitude  as 
0.00418  and  the  phase  as  -237  degrees  which  closely  agrees  with  the  Bode  plot  information 
( Table  1).  However  points  (1)  and  (2)  should  work  as  well  but  their  values  for  amplitude  and 
phase  are  0.00157  and  -211  degrees  respectively,  considerablv  different,  given  how  well  the 

other  data  set  matches  the  Bode  plots.  If  the  amplitude  and  phase  derived  from  points  (3)  and  (4) 
are  used  to  calculate  points  (1)  and  (2)  we  obtain  approximately  the  correct  values.  Obviously 
the  amplitude  and  phase  derived  from  points  (1)  and  (2)  will  not  generate  points  (3)  and  (4). 

Figure  50  is  a  plot  of  the  analytical  response  and  integration  method  output.  1  he  integrated 
output  is  about  three  times  larger  than  the  analytical  response  and  lags  it  by  50  or  62  degrees. 
The  reason  lor  the  ambiguity  in  phase  is  discussed  in  Appendix  B.  The  RK2  integration 
method  does  not  work  well  at  the  40  Hz  input  frequency  because  of  the  10  msec  integration 
step  size.  For  a  40  Hz  input  an  integration  step  of  4  msec  would  be  more  appropriate.  However 
the  digital  filters  are  still  performing  well  at  the  10  msec  sampling  interval. 

Figure  51  illustrates  the  large  spacing  of  the  sampling  intervals  for  the  40  Hz  input  frequency 
where  the  plotted  variables,  YANALS  and  YANALR,  are  the  analytical  filter  response 
calculated  at  a  10  msec  interval  and  a  2  msec  interval,  respectively.  In  Figures  49  and  50 
YANALS  were  plotted  by  drawing  straight  lines  between  data  points  rather  than  holding  a 
data  point  value  until  a  new  value  was  calculated. 

For  a  47  Hz  sine  wave  the  digital  filters  compare  as  shown  in  Figure  52.  The  sampled  47  Hz 
sine  wave  looks  like  a  modulated  signal  as  shown  in  Figure  53  where  the  invariant  and 
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Figure  48.  Plot  of  the  impulse  invariant  and  bilinear  digital  filters  with  a  10  msec 
sampling  interval. 
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Figure  49.  Plot  of  filter  input  and  filter  out,. '.3  for  the  Impulse  invariant  method  and 
the  analytical  expression  with  a  10  msec  sampling  Interval. 
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Figure  51.  Illustration  of  how  large  the  digital  sampling  interval  is  with 
respect  to  the  frequency  content  of  the  filter  output  where  the 
analytical  filter  output  was  used  for  plotting. 


Figure  52.  Plot  of  the  impulse  invariant  and  bilinear  digital  filters  -1 0  msec  sampling 
Interval. 


analyatical  responses  are  also  plotted.  Figure  54  also  contains  plots  of  the  invariant  and 
analytical  responses  but  a  bias  is  added  to  YIN  V  to  aid  in  comparing  the  two  signals.  In  Figure 
54  YINV  appears  to  be  lagging  YANALS.  However  from  Table  I  the  opposite  is  seen  to  be 
true.  Thus  the  relative  shifting  of  modulation  peaks  between  signals  is  not  an  indicator  of 
which  signal  is  leading  or  lagging  the  other.  In  general  the  digital  filters  are  doing  a  good  job 
even  though  they  are  being  updated  at  essentially  the  Shannon  limit. 

As  would  be  expected  the  RK.2  method  of  implementing  the  filter  does  not  perform  well  at 
the  10  msec  integration  step  size,  as  can  be  seen  in  Figure  55  and  Table  1.  Figure  56  illustrates 
the  sampling  coarseness  with  respect  to  the  47  Hz  frequency. 

As  previously  discussed  the  values  for  amplitude  and  phase  in  Table  l,  as  calculated  from 
the  transient  response  data  at  47  Hz,  were  based  on  peak  values. 

STOP 

BATCH,  PRINT,  PRINT,  3D,  HDMDD 

terminates  the  simulation  and  provides  a  line  printer  listing  of  the  recorded  data  for  all  of  the 
simulation  runs. 

The  Butterworth  digital  filters’  transient  response  matched  their  Bode  plot  amplitude  and 
phase  data  very  well.  The  impulse  invariant  response  closely  follows  the  continuous 
Butterworth  response  over  most  of  the  frequency  range  of  interest.  In  the  region  where  the 
exponentially  decaying  terms  have  sizeable  values  the  transient  responses  of  the  digital  filters 
do  not  overlay  the  analytical  filter  plots  but  are  indicative  of  the  analytical  response. 
Amazingly  enough  the  digital  filters  were  relatively  accurate  right  up  to  the  Shannon  limit. 

The  RK2  numerical  integration  implementation  of  the  continuous  filter  showed  an 
excellent  phase  match  and  an  adequate  amplitude  match  with  the  analytical  filter  until,  of 
course,  the  input  sine  wave  frequency  became  too  large  for  the  10  msec  integration  step  size  to 
handle.  The  RK2  algorithm  actually  outperformed  the  digital  filters  for  a  delayed  step  input-in 
the  important  sense  that  it  did  not  differ  in  phase  with  the  analytical  response. 

A  general  purpose  ACSL  macro  for  implementing  Z-transform  filters  (with  zero  I.C.’s)  was 
coded  and  tested.  Details  of  the  algorithm  are  discussed  in  Appendix  A. 
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Figure  55.  Comparison  of  the  numerical  integration  method  and  the  analytical 
response  - 10  msec  integration  step  size. 
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APPENDIX  A 

MACRO  CODE  FOR  ZXSFRM 


2/2 


AD-A140  292 
UNCLASSIFIED 


ACSL  (ADVANCED  CONTINUOUS  SIMULATION  LANGUAGE) 
SIMULATION  OF  A  THIRD  ORDE.  (U>  MITCHELL  AND  GAUTHIER 
ASSOCIATES  HUNTSVILLE  AL  E  E  MITCHELL  ET  AL.  FEB  84 
DRSMI/RD-CR-84-8  DAAK40-78-D-0010  F/G  9/1 


■■■■■■■■■■a 


NL 


! 


MACRO  MACRO  ZXSFRM (P«C) 

MACRO  REDEFINE  Y.I.YlC.JtVP 

MACRO  RELABEL  SNOOI ,SN002'SN003*SN004.SNOOS 

MACRO  ASSIGN  N 

MACRO  MULTIPLY  0 

MACRO  INCREMENT  Q(4) 

MACRO  DIVIDE  0(3) 

MACRO  IF (N*0)999 
MACRO  MULTIPLY  0 
MACRO  INCREMENT  0(4) 

ARRAY  Y (N) *YIC (N) 

CONSTANT  YIC*N*0.0 
PROCEDURAL (P ( 1 ) *P (2) *P (3) #P  (4) > 

IF (ZZFST (I ) .LT *0*5)60  TO  SN005 
DO  SNOOI  1*1 «N 

SNOOK. CALL  ZZICS(YlCfl) «¥(I>) 

YP*P (2) 

MACRO  DECREMENT  I 
DO  SN002  1*1 *N 
MACRO  INCREMENT  1 
SN002..YP*YP-P(4) (I*1)*Y(N-U 
Y(N)*YP/P(4) (1) 

YP*0.0 

MACRO  INCREMENT  1 
DO  SN003  1*1 *0(3) 

SN003. • YP*YP»P (3) (I)«Y(N-I) 

P ( 1 ) *YP 

MACRO  DECREMENT  2 
DO  SN004  1*1 *N 
SN004. «Y (I)*Y(I*1) 

SNOOS. .CONTINUE 
END 

MACRO  EXIT 

MACRO  999. .PRINT  NUMERATOR  GREATER  THAN  DENOMINATOR 
MACRO  ENO 


Figure  A.  Listing  of  the  AC8L  cods  for  macro  ZXSFRM. 
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Figure  A  is  a  listing  of  the  ZXSFRM  macro.  Its  intended  use  is  the  implementation  of  any 
order  Z-transfer  function.  The  macro  algorithm  is  based  on  the  M-method  canonical  form  of 
the  transfer  function. 

A  standard  call  to  ZXSFRM  would  appear  in  ACSL  code  as, 

ZXSFRM  (Y=X,  R,  S) 

where  y  is  the  output  of  the  filter;  X  is  the  input;  R  is  the  numerator  coefficient  array;  and  S  is 
the  denominator  coefficient  array.  The  numerator  and  denominator  coefficient  arrays  are  in 
ascending  powers  of  (Z"‘). 

Referring  to  Figure  A,  the  statement 
MACRO  MACRO  ZXSFRM(P.Q) 

obviously  does  not  correlate  with  the  argument  description  discussed  above.  So  the  following 
explanation  is  in  order.  P  is  an  array  containing  all  of  the  argument  list  variables  in  the  call  to 
the  ZXSFRM  macro.  Q  is  an  array  containing  the  first  dimensional  argument  for  each 
variable  in  the  argument  list.  For  instance,  if  R  is  an  array  of  2  and  S  is  an  array  of  3,  P  and  Q 
would  be  defined  as  follows: 

P(1)=Y 

P(2)=X 

P(3XU=R(D 

P(3X2)=R<2) 

P(4X1)=S(1) 

P(4X2)=S(2) 

P(4X3)=S(3) 

Q(l)=l 

Q(2)=l 

Q(3)=2 

Q(4)=3 

The 'MACRO  REDEFINE..' and ‘MACRO  RELABEL...’ statements  define  local  variable 
and  statement  label  names.  'MACRO  ASSIGN  N’ assigns  the  actual  number  of  arguments 
from  the  macro  call  to  the  integer  variable  N.  N  can  be  made  to  have  any  value  desired  through 
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use  of  appropriate  macro  math  operators.  N  aids  in  the  generation  of  specific  macro  code  in 
ACSL  when  an  ACSL  program  is  being  compiled. 

MACRO  MULTIPLY  0 

MACRO  INCREMENT  Q(4) 

MACRO  DIVIDE  Q(3) 

MACRO  1F(N=0)  999 
• 

MACRO  999.. PRINT  NUMERATOR  GREATER  THAN  DENOMINATOR 

is  the  code  used  to  prevent  the  attempted  implementation  of  Z  polynomials  where  the 
numerator  order  exceeds  the  denominator  order. 

MACRO  MULTIPLY  0 

multiplies  N  by  zero. 

MACRO  INCREMENT  Q(4) 

adds  the  value  of  Q(4)  to  N.  Thus  N  is  now  equal  to  the  dimension  of  the  denominator 
coefficient  array  of  the  Z-transfer  function. 

ARRAY  Y(N),  YIC(N) 

defines  the  intermediate  state  variable  array  Y  (as  for  the  M -method  formulation)  and  the 
intermediate  state  variable  initial  conditions,  YIC.  A  PROCEDU  R  AL  Mock  is  used  as  there  is 
branching  in  the  code, 

PROCEDURAL  (P(!)-P(2),P(3),P(4)). 

IF  (ZZFST(1).LT.0.5)GO  TO  SNOW 
# 

SN005..CONTINUE 

causes  the  filter  calculation  to  be  skipped  unless  this  is  the  first  minor  step  of  a  multi-step 
integration  algorithm.  ZZICS  is  called  to  initialize  the  intermediate  state  variable  array  and  to 
also  reset  the  array  to  its  last  values  when  the  REINIT  run  time  command  is  used. 


DO  SN001  1=1, N 
SN001..CALL  ZZICS(YIC(I),Y(I)) 


The  intermediate  state  variable  array  (the  (Z"1)'  M  states)  is  solved  as  shown  in  Equation  (8) 
of  Section  2, 

YP=P(Z) 

MACRO  DECREMENT  1 
DO  SN002  1=1, N 
MACRO  INCREMENT  1 
SN002..  YP=YP-P(4KI+ 1  )•  Y(N-I) 

Y(N)=YP/P(4XI) 


where 

Y(N)=M 

Y(N-i)=(Z',yM 

The  output  of  the  transfer  function  is  then  solved  by  multiplying  the  intermediate  states  by  the 
numerator  coefficients  as  shown  in  Equation  (9)  of  Section  2. 

YP-0.0 

MACRO  INCREMENT  I 
DO  SN003  1=1,  Q(3) 

SN003..YP=YP+P(3XI)*Y(N-I) 

P(1)=YP 


The  intermediate  state  array,  Y,  is  updated  in  preparation  for  the  next  calculation  interval, 

MACRO  DECREMENT  2 
DO  SN004  1=1,  N 
SN004..  Y(I)=Y(I+1) 

SN005..CONTINUE 

END 

MACRO  EXIT 

MACRO  999.. .PRINT  NUMERATOR  GREATER  THAN  DENOMINATOR 
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APPENDIX  B 

CALCULATIONS  FOR  OBTAINING  THE  TRANSIENT 
RESPONSE  AMPLITUDE  AND  PHASE  OF  THE  FOUR 
FILTER  IMPLEMENTATIONS 


From  Equation  (4)  of  Section  2,  the  steady  state  value  of  the  analytical  response  is 

D 

A  cos  7  t  +  ^  sin  7 1  (Bl> 

or 

P  (sin  <j>  cos  7  t  +  cos  <|>  sin  7t) 

*  P  sin  ( 7t  +  (j» 

where 

F=the  steady  state  amplitude  of  the  filter  output 
^=phase  shift  of  the  filter  output  with  respect  to  the  input 
<fy<0  is  phase  lag 
4>X)  is  phase  lead 
So 

<j)  -  tan-1  (~) 

where  A  and  B  are  defined  in  Equation  (4)  of  Section  2  and 

P  ■  A/sin  $  or  B/(7cos((»  (B2) 

Due  to  the  arbitrary  quadrant  location  of  the  arctangent  function  as  shown  below, 

-  |  +  tan 

tan  (8  +  mr)  ■  same  value, 

n  »  0,  1,  2,  ... 


8 


the  phase  lag  from  the  Bode  plots  and  a  corresponding  rough  measure  of  phase  lag  in  the 

%4 

transient  response  plots  are  used  as  justification  for  the  choice4^  made.  There  are  also  certain 
restrictions  on  4».  Substituting  in  ( I )  for  A  and  B  the  actual  terms. 


<j>  =  tan 


(y  <&2  -  2  \ 


(B-3) 


which  except  lor  the  5  Hz=y  case  <7=10.  20.  40....  11/ 1  is  always  going  to  be  less  than  zero. 
Secondly  the  filter  will  have  phase  lag  and  therefore  <£<0.  Combining  these  two  facts  with  the 
knowledge  that  |d>|<360  degrees  leaves  only  two  choices  for  <f>:  tan  '(...)  or  tan  1  (...)-7r. 

■  •  me  true  analytical  values  for  amplitude  and  phase  can  be  determined.  The  Bode  plot 
- 1 .  .i  three  ot  the  four  filters  may  be  obtained  from  a  modified  version  of  the  program  in 
•».  ,  non  3.  I  he  amplitude  and  phase  values  for  the  two  digital  filters  and  the  integration 
method  can  be  obtained  from  the  transient  response  program  data. 


1  or  a  linear  transfer  function  the  steady  response  to  a  sine  wave  is  another  sine  wave  shifted 
in  phase  and  changed  in  amplitude, 

sin  tat  =>  A  sin  (<ot  +  <j>)  (B-3) 


where,  as  noted  above,  dKO  is  phase  lag  and  <£X)  is  phase  lead.  For  two  arbitrary  points  Xi, 
and  X:,  in  the  transient  response  output  we  have 


A  sin(u>t^  +  <J>)  * 

A  8in(ut2  +  4>)  ■  X2 
where  from  (4) 

-i  xi 

<t>  -  sin  1  (£=■)  - 


(B-4) 

(B-5) 


(B-b) 


Substituting  (6)  into  (3)  gives 


Letting  >=cot2-<oti, 


Ajsln  y  cos^sin”1  -—-^+  cos 


or 


Va2-  xj 

A  I  sin  y  - s - +  cos  Y 


I 


ft-)] 


X2. 


Squaring  both  sides  and  solving  for  A, 


’  /  X2  -  X ^cos  Y  \ 

2  2l 

\  sin  y  ) 

(B-7) 


where  the  +  sign  is  chosen.  As  mentioned  above  there  is  a  multitude  of  angles  which  satisfy  (6). 
Due  to  the  arbitrary  quadrant  location  of  the  arc  sine  function  as  shown  below. 


0  +  2wr 

SIGN  (180,  0)  +  2nir  -  6 
n  *  0,  +  1,  +  2,  ... 


the  phase  lag  from  the  Bode  plots  and  a  corresponding  rough  measure  of  phase  lag  in  the 
transient  response  are  used  to  aid  in  determining  4>.  Again  -360<^p  so  if  4>  is  calculated 
positive  -2ir  is  added  to  it  and  if  4<0  it  can  only  have  two  values:  4>  or  -180-|<£. 

Equation  (B-7)  is  undefined  when  y,  or  wti-wti,  is  a  multiple  of  it.  Other  quirks  in  (7)  were 
found  and  are  mentioned  for  the  appropriate  cases  in  Section  4.  But  evidently  the  problem 
stems  from  the  fact  that  cut2-wti,  can  also  be  interpreted  as  another  w*t2*-«u*ti*  as  long  as  the 
value  of  the  subtraction  remains  the  same.  However  the  amplitude  expression  holds  quite  well 
at  the  lower  frequencies  and  for  the  purposes  of  this  report  were  considered  adequate.  ti,  t2,  Xt 
and  X2  are  judiciously  picked  near  the  end  of  the  plotted  filter  responses  to  obtain 
approximately  steady  state  values. 
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